#include "dEdxGausN.h"

Double_t fres = -999;

Double_t Fproton(){   return 0.20;}
Double_t Fpion(){     return 0.49;}
Double_t Fkaon(){     return 0.30;}
Double_t Felectron(){ return 0.01;}

Double_t TotalN(const Double_t xx)
{
  const Double_t ex = -5/3. * xx + 13./3. -1;
  return pow(10, ex);
}

TRandom3 fran(123);

Double_t bbpar[10];
//const Double_t mBBIni[]={  3.43178e+01, 2.40223e+00,  1.90810e-01, 2.70186e-01, 7.40610e+01};
//new para with p2 : K-pi crossing momentum in GeV/c, p3: K-p crossing momentum in GeV/c
const Double_t mBBIni[]={3.43178e+01, 2.40223e+00,  1.06824287379018523e+00, 2.54489579623403639e+00,  7.40610e+01};
/*
gSystem->Load("lib/libBaseUtils.so");
gSystem->Load("lib/libTPCHist.so");
gSystem->Load("lib/libdEdxGausN.so");
const Double_t tmp[]={3.43178e+01, 2.40223e+00,  1.90810e-01, 2.70186e-01, 7.40610e+01};  
dEdxGausN::FindBetheBlochCrossing(tmp, dEdxGausN::Mkaon(), dEdxGausN::Mpion())
dEdxGausN::FindBetheBlochCrossing(tmp, dEdxGausN::Mkaon(), dEdxGausN::Mproton())
const Double_t tmp2[]={3.43178e+01, 2.40223e+00,  1.06824287379018523e+00, 2.54489579623403639e+00,  7.40610e+01};
TF1 *ff=new TF1("ff",dEdxGausN::myBetheBloch, 0.1, 1e5,5)
ff->SetParameters(tmp2)
ff->SetLineColor(kBlue); ff->Draw()
ff->SetMaximum(200); gPad->SetLogx()
ff->SetMinimum(0); ff->SetMaximum(200); gPad->SetLogx()
TF1 *f2=new TF1("f2",dEdxGausN::myOldBetheBloch, 0.1, 1e5,5)
f2->SetParameters(tmp)
f2->Draw("same")
 */
const Double_t alephIni[]={1.63568e+00, 2.20978e+01, -2.43255e+01, 2.03191e+00, 3.97478e+00};
const Double_t lundIni[]={ 3.36603e+01, 1.97986e-01,  2.34633e+01, 9.91833e-01, 7.40744e+01};

const Double_t fpmin = 0.15;
const Double_t fymax = 150;

Double_t FillHist(TH2D * hh, TH2D *mcmom, const Int_t ibin, const Double_t frac, const Int_t itype)
{
  if(fres<0){
    printf("FillHist wrong fres %e\n", fres);exit(1);
  }

  const Double_t mass = dEdxGausN::Mass(itype);

  const Double_t xx = hh->GetBinCenter(ibin);
  const Double_t mom = pow(10, xx);
  if(mom<fpmin){
    return 0;
  }

  const Double_t nsample = TMath::Max(200.0, TotalN(xx))*frac;

  //printf("mass %e nsample %.0f totaln %e frac %e\n", mass,nsample,TotalN(xx), frac);;

  const Double_t my = dEdxGausN::GetMean(xx, mass, bbpar);
  //printf("ibin %d xx %e mass %e my %e bbpar %e %e %e %e %e\n", ibin, xx, mass, my, bbpar[0], bbpar[1], bbpar[2], bbpar[3], bbpar[4]);

  for(Int_t ip=0; ip<nsample; ip++){
    const Double_t val = fran.Gaus(my, my*fres);
    if(val>fymax || val<0)
      continue;
    //printf("fill xx %e val %e my %e fres %e\n", xx, val, my, fres);
    hh->Fill(xx,val);
    mcmom->Fill(xx, itype);
  }

  return nsample;
}

void genhist()
{
  TList *ll=new TList;

  TH2D * horiginal = new TH2D("horiginal","ToyMC", 300, -1, 2, 1200, 0, 200);
  ll->Add(horiginal);

  TH2D *mcmom = new TH2D("mcmom", "trueMC", 300, -1, 2, 5, -0.5, 4.5);
  ll->Add(mcmom);

  if(dEdxGausN::BBOpt()== dEdxGausN::kmBB){
    memcpy(bbpar, mBBIni, sizeof(Double_t)* dEdxGausN::NparBB());
  }
  else if(dEdxGausN::BBOpt()== dEdxGausN::kRawALEPH){
    memcpy(bbpar, alephIni, sizeof(Double_t)* dEdxGausN::NparBB());
  }
  else if(dEdxGausN::BBOpt()== dEdxGausN::kLund){
    memcpy(bbpar, lundIni, sizeof(Double_t)* dEdxGausN::NparBB());
  }
  else{
    printf("FillHist BBOpt wrong ! %d\n", dEdxGausN::BBOpt()); exit(1);
  }

  for(Int_t ii=0; ii< dEdxGausN::NparBB(); ii++){
    printf("inputpar %d %e\n", ii, bbpar[ii]);
  }

  Double_t ntot=0;
  for(Int_t ibin = 1; ibin<=horiginal->GetNbinsX(); ibin++){
    ntot += FillHist(horiginal, mcmom, ibin, Fproton(),   dEdxGausN::kPROTON);
    ntot += FillHist(horiginal, mcmom, ibin, Fpion(),     dEdxGausN::kPION);
    ntot += FillHist(horiginal, mcmom, ibin, Fkaon(),     dEdxGausN::kKAON);
    ntot += FillHist(horiginal, mcmom, ibin, Felectron(), dEdxGausN::kELECTRON);
  }

  printf("ntot %e\n", ntot);

  TFile *fout = new TFile(Form("outgenhist_TOY%d.root", dEdxGausN::Ntype()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=2){
    printf("argc!=1\n");return 1;
  }

  fres = atof(argv[1]);
  printf("fres = %e;\n", fres);

  dEdxGausN::Ini();

  genhist();
  return 0;
}
